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We use the one-dimensional steady version of the equations derived in paper I to compute the 
, structure of shock waves. The agreement with experiment is good, especially when we retain the 

' experimental value of the Prandtl number adopted in II. 

■ I. INTRODUCTION 
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In the simplest descriptions of shock waves, one uses the Euler equations of fluid dynamics to develop jump conditions 
across the discontinuity that describes the shock. These conditions have qualitative value but there are problems for 
which such limiting solutions are not adequate. For example, in studying radiation from shock waves, say to infer the 
: properties of the radiating atoms, the conditions within the shocks may become important. If the matter ahead of 
J> the shock is neutral and that behind the shock is fully ionized, the observed spectral lines may be formed in the shock 
OO itself and so conditions there need to be carefully worked out. This example is one of many that we might have cited 
\l to motivate our present study of the structure of shock waves. Even a cursory look at the literature on this problem 
gives a clear picture of the importance that has been attached to it. 

Much of the work on the structure of shock waves has been aimed at improving on the description provided by the 
Navier-Stokes equations. A review of such attempts was given by Galkin et. al. M, who compared various higher 
order solutions of Boltzmann equation with results from the Chapman-Enskog Q and Hilbert methods jlj] . In 
terms of shock structure alone, the higher order solutions give a significantly improved result |Q, but a number of 
fundamental problems concerning the status of these equations remain open H, ^ It is especially troubling 
that a large number of higher order nonlinear terms in some of the proposed improvements make it difficult to use 
O I the results in realistic problems. Similar difficulties beset extensions of Grad's moment method ||l3|, such as 
the Extended Irreversible Thermodynamics p^ , which does give a good result for the problem of ultrasonic sound 
^ ' wave propagation. It remains true however, that a large number of equations must be solved to achieve a reasonable 
accuracy with EIT. 

When the densities arc sufficiently low, direct numerical simulation by the Monte Carlo method is the most reliable 
way to compute high-Knudsen-number flows, though the computational cost may be high in regimes near continuum 
limits A relatively new and effective method that has been described for direct solution of the relaxation (or 

cd BGKW) model does offer hope of improvement but it remains to taken beyond the two-dimensional problems 
that it has so far tested very well on. 

It thus appears that an effective macroscopic description would be of value for such problems and so, in this paper, 
we show how the equations derived in paper I may fill this need. To do this, we work out the structure of shock 
waves in one-dimensional steady flow starting from the equations of paper I. Those equations were derived from 
kinetic theory without using some of the traditional simplifications associated with the Chapman-Enskog approach. 
In particular, we did not use results from lower order approximations to simplify higher order equations. In paper II, 
the first-order development (in mean free path) of paper I was tested against observations of ultrasound propagation. 
We found good quantitative agreement between our theory and the experiments in the high Knudsen number regime 
where Navier-Stokes equations clearly fail. Though the success of our approach works well beyond the expected limits 
of validity of the theory, such occurrences are not unheard of in good asymptotic methods. 

In the present paper, we go on to see how well the theory works for the computation of the structure of shock waves. 
Since the shock thickness is typically of the order of a mean free path, we again are pushing against the limits of 
validity of the expansions used in deriving the fiuid equations. We confront in addition the added challenge of strong 
nonlinearity. Moerover the fiow considered in shock theory can is far from thermodynamic equilibrium, so this too 
makes for a stringent test. 

A simplification that makes the comparison relatively easy to draw is that, in the study of shocks, we may separate 
the continuum differential equations of fiuid motion from the boundary conditions that must be stated to complete 
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a well-posed problem. The boundary conditions for the shock problem is not in question since we may impose the 
aforementioned jump conditions, expressed in terms of the Rankine-Hugoniot relations, to serve in the role of the 
boundary conditions. Hence in this study of shock structure, we are able to focus on effects due only to the differential 
equations themselves and so examine the validity of our version of the fluid equations. 

In Sec. II, we recall the fluid dynamical equations derived in paper I. Then, in Sec. Ill, the structure of one- 
dimensional shock waves are computed and compared to analogous results obtained with the Navier-Stokes equations 
as well as with experiments. The paper concluded with a brief discussion in Sec. V. 



II. STATEMENT OF THE EQUATIONS 



In paper I we proposed a modification of the usual asymptotic techniques for deriving fluid equations from kinetic 
theory. Our procedure avoids the simplification introduced by Chapman and Enskog in which the results of lower 
approximations are introduced into the higher approximations. When we proceed in this way, we obtain these fluid 
equations ||, |6| : 



dtp + V ■ (pu) = 
dt u + \i-\7\x+-\7 



(2.1) 
(2.2) 



dtT + u-VT 



3pR 



(P: Vu + V-Q) = 0, 



(2.3) 



where p is the mass density, u is the average velocity of the particles in a fluid element cell, and T is the temperature. 
We assume that the gas is made up of identical, structureless particles with mass m so that R = k/m is the gas 
constant with k the Boltzmann constant. Our expression for the stress tensor is 



where 



and 
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pE : VVu 



p = TpRT 

is the viscosity expressed in terms of the mean free time r. For the stress tensor we found 

(9u* 



and for the heat current, we have 



Q = -ryV ln(pr" 



— — TiVlnT Uir— 

2' 2^Dt 



where rj — ^pRT for the relaxation model. 

By contrast, in the Navier-Stokes equations with no bulk viscosity, one has 

P = pRTI - pE : VVu 

and 



(2.4) 



(2.5) 



(2.6) 



(2.7) 



(2.8) 



(2.9) 



Q = -T]V InT 



(2.10) 
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III. SHOCK THEORY 



A. Basic equations 



As inidcated in Fig. y| when the velocity of the flow in the upstream (a: — oo) exceeds the sound speed of the 
medium, a shock front forms. In a frame comoving with the shock front, we see a steady shock layer, whose structure 
is determined by the upstream thermodynamic quantities and the dissipation mechanism. The structure of the shock 
wave provides a straightforward test of the equations. 



shock layer 
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FIG. 1: The schematic diagram of a shock wave 



Let us use subscripts 1 and 2 to denote values at large distances upstream (x — > — oo) and downstream (x — *■ oo) 
of the shock front, with u = {v{x) , 0, 0), p and T as functions of only x. In the frame comoving with the shock front, 
we have dt ~ 0- On integrating to (2.3) from uniform upstream state to an arbitrary position x in the shock, 
we obtain 



pv = piVi 
p + pv^+E^^ 



Pi + PlVl 



pv{CyT + -W^) +pv + VTTxx + Qx = PlVi{CyTi + -vl) +piVi 



(3.1) 
(3.2) 

(3.3) 



where c„ = En is the xx component of the viscous stress tensor, and is the x component of the heat current. 
The simple form of the right sides of these equations results from the vanishing of the derivatives of the fluid variables 
far from the shock. 



Prom (3.1), we can express p in terms of ^ 



PlVl 



(3.4) 



Combination of ( p.4| ) with (3_^) and ( |3.3| ) leads to 



= Rpivi ( — I + pivi{vi - v) 

\Vl V J 

3 1 
= -RpiVi{Ti - T) + pi{vi - v) + -pivi{vi 



(3.5) 
(3.6) 



We nondimensionlize these equations using vi as the unit of speed Ti as unit of temperature and the mean free 
path Ai at upstream infinity as the unit of length. Then, we introduce the nondimensional quantities X = x/Xi, 
w{X) =v/vi, e{X) = T/Ti and 



w{X) 



E^ 

Pi 



and 



PlVl 



(3.7) 
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The equations become 



where the Mach number is 



and 7 = 5/3. 



w — 1 — 
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(3.8) 
(3.9) 
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with 
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We see in from (2.6) that the viscosity is /i = rp, where r is the mean flight time of the particles. In turn, r is 
the mean free path over the mean speed, which is the local speed of sound. Here, we adopt the simplest form of the 
relaxation model, namely that with constant r. Then, we follow Gilbarg et al. and take 



M = 



(3.10) 



where s is a parameter that is adjusted according to the nature of the constituent particles. For argon, s = 0.816 
has often been used ||], but Kestin et al. 17| suggest that s — 0.64 is a better value for argon. We shall adopt the 
latter, more recent value in these computations. As to the conductivity, the relaxation model gives r] = ^^RT but, 
when the Boltzmann collision term is used, we obtain a slightly different value in better agreement with experiment. 
As the difference results from the atomic model rather than to fluid dynamical issues, we shall adopt the formula 
?7 = 15/ii?r/4 to remove the effects of the inaccuracy of the atomic model. 

To these formulae we adjoin the the closure relations (2.4) and (2.S ) that may be rewritten in nondimensional form 

as 



vu = — A / - Ml s 



3 5s 
* ^ ^2V6Mi 
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29' + —Mlww' - e— 
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(3.11) 
(3.12) 



where the primes indicate differentiation with respect to X. 

On combining (3.11) and (3.12) with (3.8) and (3.9), we obtain these equations: 
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(3.13) 
(3.14) 



For comparison, we also write down the equations determining the shock structure from the Navier-Stokes equations: 
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4 /5 
sV 6 
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-+ Ml(l-w) 
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+ ^-ii-w) + lM!{i-wr 



(3.15) 
(3.16) 



B. Critical Mach number 



Before carrying out the numerical integration of equations ( 3.13| )— (3.14), we examine them to draw some general 
conclusions. First, we o bserv e that, since the x-der ivati ve s of it; and 9 vanish far upstream and downstream, the 
right sides of equations ( ^.13 )-(3.14) (as well as of (3.15)-( p.l6[ )) must vanish there. These conditions provide two 
simultaneous equations which are readily solved for the fixed points {wi, 9i) =^ (1,1) and 



(^2,6*2) = 



Ml 



iMf 



3 5Mf + UMl - 3 
16Mf 
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Our aim is to find solutions that connect the two fixed points. However, in the case of our equations, though the 
fixed points must occur at the locations just found in the {w,9) plane, those locations need not be fixed points: the 



left sides of (3.13)-(3.14) may vanish when the determinant of the matrix of the coefficients of the derivatives vanishes. 



This is seen when we write the equations in the standard form 

w'\ _ 1 [6 fvj/M_ 



where 



with ^ = lOMf/9. 

If we can solve (3.17)-( pJ8 ) to obtain explicit expressions for the derivatives, we can then solve the coupled first- 
order ordinary differential equations to find solutions connecting the fixed points. This is possible if the determinant 
of M does not vanish. The critical condition is then obtained by setting the determinant to zero, a step that leads to 
the relation 

= ^M>2 . (3.19) 
9 

At upstream infinity, we have {w, 9) — (1, 1) and this leads to the critical value Mi — 3/V2 « 2.12. We cannot solve 
for the derivatives when Mi exceeds this value and so the theory breaks down above this value. 

When the upstream flow speed is large enough to make the Mach number surpass the critical Mach number Mf , 
D changes its sign. If this leads to either > or < near the upstream, which is true for our case here, 
it will be impossible to match the downstream values monotonically. This happens for any hyperbolic system, as 
is commonly seen in moment formalisms [|l^ and first noticed in Grad's 13-moment method. Interestingly, it is 
numerically verified [|l6|, |l^ that only when the largest upstream critical Mach number is surpassed, does the above 
situation happen. As more and more moments are included, the largest critical Mach number becomes larger and 
larger, which in a sense partly relaxes the constraints on the application to shock study with moment method. It 
is thus expected that with the higher-order terms included in our formalism, the range of the Mach number will be 
extended. Nevertheless, the constraints on the Mach number is not due to the expansion scheme, but due to the 
nature of the relaxational model. In a following paper, we shall show that, when classical Boltzmann integral is used 
in place of the relaxation term, the constraints on Mach number will be removed. However, physically, information 
speed can not be infinite, which indicates that the final version of hydrodynamics must be of hyperbolic nature. 
Intuitively, this issue can only be addressed consistently under the frame of relativity. Amazingly, when our modified 
Hilbert expansion is applied to relativistic Boltzmann equation a causal hydrodynamics which has constraints 
on the Mach number originated from relativity is obtained. Such a hydrodynamics is free of the stability problem 
and causality problem inherent in the Steward-Israel version of hydrodynamics |l^, ^ , the relativistic counterpart 
of N-S formalism. For our present purpose, we are content with the nonrelativistic version. The major concern here 
is whether the accuracy is improved when Knudsen number is not small. 

C. Comparisons 



We integrate the equations ( 3.13 )^( 3.14 ) with different Mach number and make comparisons on the following 
aspects. 

1. Shock profiles 

In order to facilitate the comparison, let us first introduce quantities dP and dQ which are defined as the difference 
of stress and the heat current between the new version and the conventional N-S formalism. The definitions are the 
following: 



Vo InT-f ^V-u 



I , (3.20) 
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For one-D flow, we obtain 



with AP and AQ given by 



AQ = - (?7Vln(pT) + -/iPou 



AP = APxx , 
AQ = AQx 



AP = -^i{VolnT+-v'] , 



AO = - |7?[ln(pr)]' + ^Pof I 



(3.21) 

(3.22) 
(3.23) 

(3.24) 
(3.25) 



In the study of shock wave, the upstream quantities are used to rewrite the aforementioned dimensionless dP and 



dQ: 



dP: 



AP 
Pi 



sMi 



-cu 



(3.26) 



dQ: 



2AQ 

3pivi 



sMi 



156 Lu 



(3.27) 



From the above definitions, we see that dP is equivalent to the bulk viscous pressure and dQ accounts for the extra 
heat current due to the non-adiabatic effect, both of which are taken as zero in the derivation of N-S formalism. Such 
an approximation is not justified when Knudsen number is not small, as will be shown in the following. 
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FIG. 2: Density profiles with Mach number M=1.5 

Figures show the density, velocity and temperature profiles with Mi — 1.5. It is seen that the profiles with the 
modified hydrodynamics are wider than those with N-S formalism. Such a trend is a well-known fact and has been 
verified experimentally. 

The difference in the profiles are due to the nonzero dP and dQ, plotted in Fig. ^ and Fig. ^, respectively. A simple 
estimation shows that the dissipation due to these two extra effects increases monotonically from the downstream to 
the upstream, which tends to lower the gradients of the thermodynamical variables and results in a wider shock. Such 
effects become more noticeable when Mach number is larger, as is shown in Figs. [Tj^'-^pl], where the same quantities as 
those in Figs. ||~^ are plotted but with Mi = 2. Compared to the case with Mi — 1.5, the profiles are steeper and 
more asymmetric. The difference of the two versions of hydrodynamics is now more significant due to the increase in 
the Knudsen number with the Mach number. From Figs. p^|^ [Tl| , it is clear that the bulk viscous pressure and the 
non-adiabatic heat current are no longer negligible. 
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FIG. 3: Velocity profiles with Mach number M=1.5 




FIG. 4: Temperature profiles with Mach number A/=1.5 



2. Inverse thickness, asymmetry and T-p separation 



The inverse thickness is defined as 



5-' 



P2 - Pi 



(3.28) 



where p'^ is the maximum value of the gradient of density profile, and as before pi and p2 are the upstream and 
downstream values, respectively. 

The asymmetry parameter B is defined as 



As - 



(3.29) 



with Pa = \{p2 - Pi)- 

Finally, the temperature-density separation ApT is defined as the distance between the medium value points of 
temperature and density. 

Figure ^ shows the inverse thickness versus Mach number. The experimental data for Argon and Xeron are also 
shown. Concluded from the plot, the modified hydrodynamics results in a wider shock for a given Mach number. In 
comparison with the experimental data and the direct Monte Carlo simulation results(for full Boltzmann equation). 
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FIG. 5: dP profiles with Macli number M=1.5 
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FIG. 6; dQ profiles with Mach number 

the accuracy is significantly improved. The physical explanation for the wider shock is as follows. The process- 
dependent bulk viscosity causes extra dissipation, which, combined with the non-adiabatic effect, makes the gradients 
of thermodynamic variables smoother and thus widens the shock. 

As for the asymmetry parameter, in the case of Mi = 2.0, N-S formalism gives Ag = 1.25, the modified one gives 
As = 0.83, which is closer to the experimental data 0.93. Aside from the quantitative improvement, the modified 
version also becomes qualitatively consistent with the experiment, i.e. As should be smaller than unity. 

For the same Mach number Mi = 2.0, N-S leads to a temperature-density separation of 1.04, whereas the modified 
version leads to 1.67 and the experimental data for Argon is 1.50 

IV. CONCLUSIONS AND DISCUSSIONS 

We have showed that the modified Hilbert expansion leads to a hydrodynamics of higher accuracy than that from the 
classic Chapman-Enskog expansion. In the study of rapid variation such as high-frequency sound wave propagation, 
this modified hydrodynamics provides a result consistent with experiment. The fitting is uniformly good for all range 
of Knudsen number. For the study of shock wave structures, within the range of its critical Mach number, it provides 
results comparable to those from direct Monte Carlo simulation. The physics underlying the new expansion scheme is 
that it relaxes the requirement that the gradient terms be small as is the case in CEE. This makes the method capable 
of deriving macroscopic equations valid for processes far from equilibrium. As we have seen in the shock structure 
study, the extra dissipation caused by the new mechanism of gradient of density smoothes the density profile. Such 
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FIG. 7: Density profiles with Mach number M=2.0 
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FIG. 8: Velocity profiles with Mach number M=2.0 

a non-adiabatic dissipation is omitted in the Navier-Stokes formalism due to the extra constraints enforced by the 
CEE procedure in the derivation. As is discussed in section III, the existence of the critical Mach number is caused 
by the relaxation model we adopted for the kinetic equation. In one of our following papers [ p6[ we will show that, 
when the Boltzmann collisional integral is used in place of the relaxation term, such a constraint on Mach number 
is removed. However, a consistent hydrodynamics satisfying causality principle yet without extra constraint beyond 
relativity on the possible Mach number can only be formulated with the relativistic Boltzmann equation, as will be 
shown in Ref. [ p7| . 
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